function [DivorceCostsM,DivorceCostsF,DivorceCostsMF,solved] = get_indices(couple,singlemale,singlefemale,leisurem,leisuref,hworkm,hworkf,wagem,wagef,q,Q,nonlabor,mutually_acceptable,typem,typef,wagem_bd,wagef_bd,wbound,text_solver)

% ================
% define constants
% ================

T = 112;                                              % total potential working hours
H = length(leisurem);                                 % total number of households in the data
alpha = 0.4;                                          % bound for non-labor income

% =========================
% define decision variables
% =========================
yalmip('clear')

qM = sdpvar(H,1);                                     % private consumption of Hicksian good male
qF = sdpvar(H,1);                                     % private consumption of Hicksian good female

pM = sdpvar(H,H,'full');                              % personalized public good price male
pF = sdpvar(H,H,'full');                              % personalized public good price female

pmhworkM = sdpvar(H,H,'full');                        % personalized male's hwork price male
pmhworkF = sdpvar(H,H,'full');                        % personalized male's hwork price female

pfhworkM = sdpvar(H,H,'full');                        % personalized female's hwork price male
pfhworkF = sdpvar(H,H,'full');                        % personalized female's hwork price female

wageM = sdpvar(H,1);                                  % shadow wage male
wageF = sdpvar(H,1);                                  % shadow wage female

nl_incomeM = sdpvar(H,1);                             % non-labor income male
nl_incomeF = sdpvar(H,1);                             % non-labor income female

divorce_costsM = sdpvar(H,1);                         % stability indices for individual rationality male
divorce_costsF = sdpvar(H,1);                         % stability indices for individual rationality female
divorce_costsMF = sdpvar(H,H,'full');                 % stability indices for no blocking pairs


%%
% ==================
% define constraints
% ==================

Constraints = [];

% restriction on individual private consumption
for i=1:H
    Constraints = [Constraints, q(i) == qM(i) + qF(i), qM(i) >=0, qF(i) >= 0];
    
    if singlefemale(i) == 1
        Constraints = [Constraints,qM(i) == 0];
    elseif singlemale(i) == 1
        Constraints = [Constraints,qF(i) == 0];
    end
    
end



% shadow price should be non-negative and equal to observed wage if the individual is employed
for i=1:H    
    
    Constraints = [Constraints,wageF(i) >= 0, wageM(i) >= 0];
    if ~isnan(wagem(i))
        Constraints = [Constraints, wageM(i) == wagem(i)];
    else
        if typem(i) ~= 0
            x = find(wagem_bd(:,3)== typem(i));
            Constraints = [Constraints, wagem_bd(x,1)- wbound*wagem_bd(x,2) <= wageM(i), wageM(i) <= wagem_bd(x,1) + wbound*wagem_bd(x,2)];
        end
    end
    
    
    if ~isnan(wagef(i))
        Constraints = [Constraints, wageF(i) == wagef(i)];
    else
        if typef(i) ~= 0
            x = find(wagef_bd(:,3)== typef(i));            
            Constraints = [Constraints, wagef_bd(x,1)- wbound*wagef_bd(x,2) <= wageF(i), wageF(i) <= wagef_bd(x,1) + wbound*wagef_bd(x,2)];
        end
    end
end



% restriction on public prices
for i=1:H
    for j=1:H
        
        Constraints = [Constraints, 1 == pM(i,j) + pF(i,j), pM(i,j) >=0, pF(i,j) >= 0];
        Constraints = [Constraints, wageM(i) == pmhworkM(i,j) + pmhworkF(i,j), pmhworkM(i,j) >=0, pmhworkF(i,j) >= 0];
        Constraints = [Constraints, wageF(j) == pfhworkM(i,j) + pfhworkF(i,j), pfhworkM(i,j) >=0, pfhworkF(i,j) >= 0];
        
        if (couple(i) == 1 || singlemale(i) == 1) && singlemale(j) == 1
            Constraints = [Constraints, 1 == pM(i,j), wageM(i) == pmhworkM(i,j), wageF(j) == pfhworkM(i,j)];
        end
        
        if (couple(j) == 1 || singlefemale(j) == 1) && singlefemale(i) == 1
            Constraints = [Constraints, 1 == pF(i,j), wageM(i) == pmhworkF(i,j), wageF(j) == pfhworkF(i,j) ];
        end
    end
    
end


% feasibility restrictions on the individual nonlabor income (for each matched pair)
% individual nonlabor incomes after divorce must lie between 40% and 60% of total nonlabor income under marriage
for i=1:H
    
    Constraints = [Constraints, nonlabor(i) == nl_incomeM(i) + nl_incomeF(i)];
    
    if couple(i) == 1
        if nonlabor(i)>=0
            Constraints = [Constraints, alpha*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= (1-alpha)*nonlabor(i), alpha*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= (1-alpha)*nonlabor(i)];
        else
            Constraints = [Constraints,(1-alpha)*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= alpha*nonlabor(i), (1-alpha)*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= alpha*nonlabor(i)];
        end
    elseif singlefemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeF(i)];
    elseif singlemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)];
    end
    
end


% individual rationality
for i = 1:H
    if couple(i) == 1
        Constraints = [Constraints, wageM(i)*T + nl_incomeM(i) - divorce_costsM(i) <= qM(i) + wageM(i)*leisurem(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];
        Constraints = [Constraints, wageF(i)*T + nl_incomeF(i) - divorce_costsF(i) <= qF(i) + wageF(i)*leisuref(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];
    end
end


% no blocking pairs
for i = 1:H
    for j = 1:H
        if mutually_acceptable(i,j) == 1
            Constraints = [Constraints, (wageM(i) + wageF(j))*T + nl_incomeM(i) + nl_incomeF(j) - divorce_costsMF(i,j) <= ...
                qM(i) + qF(j) + wageM(i)*leisurem(i) + wageF(j)*leisuref(j) + pM(i,j)*Q(i) + pF(i,j)*Q(j) + ...
                pmhworkM(i,j)*hworkm(i) + pmhworkF(i,j)*hworkm(j) + pfhworkM(i,j)*hworkf(i) + pfhworkF(i,j)*hworkf(j)];
        end
        
    end
end



% bound divorce cost
for i=1:H
    Constraints = [Constraints,0 <= divorce_costsM(i), 0 <= divorce_costsF(i)];
    
    for j=1:H
        Constraints = [Constraints, 0 <= divorce_costsMF(i,j)];
    end
end


%%
% =================
% Solve the problem
% =================

% Specify the objective function
Objective = sum(divorce_costsM) + sum(divorce_costsF) + sum(sum(divorce_costsMF));
options = sdpsettings('verbose',0,'solver',text_solver);

% Solve the problem
% use -Objective if we want to maximize.
sol = optimize(Constraints,Objective,options);

solved = sol.problem;
% Analyze error flags
if sol.problem ~= 0
    sol.info;
    yalmiperror(sol.problem);
end

%get the value out
DivorceCostsM = value(divorce_costsM);
DivorceCostsF = value(divorce_costsF);
DivorceCostsMF = value(divorce_costsMF);

return
